Europhysics Letters 



PREPRINT 



A numerical renormalization group study of laser induced 
freezing 

Debasish Chaudhuri(*) and Surajit Sengupta(**) 

Satyendra Nath Bose National Centre for Basic Sciences - Block-JD, Sector-Ill, Salt 
Lake, Calcutta - 700098. 



PACS. 64.70.Dv - Solid-liquid transitions. 

PACS. 64.60.Ak - Renormalization-group, fractal, and percolation studies of phase transi- 
tions. 

PACS. 82.70.Dd - Colloids. 



Abstract. - We study the phenomenon of laser induced freezing, within a numerical renormal- 
ization scheme which allows explicit comparison with a recent defect mediated melting theory. 
Precise values for the 'bare' dislocation fugacities and elastic moduli of the 2-d hard disk system 
are obtained from a constrained Monte Carlo simulation sampling only configurations without 
dislocations. These are used as inputs to appropriate renormalization flow equations to obtain 
the equilibrium phase diagram which shows excellent agreement with earlier simulation results. 
We show that the flow equations need to be correct at least up to third order in defect fugacity 
to reproduce meaningful results. 



Introduction. ~ Re-entrant "laser-induced" freezing [1,2] (RLIF) has received consider- 
able attention [3-13] in recent times. A static interference pattern obtained by two crossed 
laser beams provides an external potential periodic in one dimension (1-d) which induces a sys- 
tem of dielectric colloidal particles, confined in two dimensions (2-d), to freeze. Surprisingly, a 
further increase in potential strength causes a reentrant melting [2] transition. Qualitatively, 
starting from a liquid phase, the external periodic potential immediately induces a density 
modulation, reducing fluctuations which leads to solidification. Further increase in potential 
confines the system to decoupled 1-d strips. The dimensional reduction now increases fluc- 
tuations remelting the system. The early mean field theories, namely, Landau theory [1] and 
density functional theory [3] predicted a change from a first order to continuous transition 
with increase in potential strength and failed to describe the reentrant behavior, a conclu- 
sion seemingly confirmed by early simulations [6]. Quite generally however, in 2-d, the effect 
of thermal fluctuations is substantial and mean field theory may obtain qualitatively wrong 
results especially for fluctuation driven transitions. Following the defect mediated disorder- 
ing approach of Kosterlitz and Thouless [14] (KT), Frey, Nelson and Radzihovsky [5](FNR) 
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Fig. 1 - The phase diagram of the hard disk system in the presence of a 1-d, commensurate, periodic 
potential in the packing fraction (n) - potential strength (/3Vo) plane. The lines in the figure are a 
guide to the eye. The dashed line denotes earlier Monte Carlo simulation results[7] and the solid line 
is calculated through our numerical renormalization group study. The dash-dotted line at rj ~ .705 
denotes the calculated asymptotic phase transition point at f3Vo = oo. 



proposed a detailed theory for the reentrant transition based on the unbinding of disloca- 
tions with Burger's vector parallel to the line of potential minima. More accurate simulation 
studies [7-13] on systems of hard disks [7], soft disks [8], DLVO [9,10] etc. confirmed the 
reentrant freezing-melting transition in agreement with experiments [2] and FNR theory. A 
systematic finite size scaling analysis [7] of simulation results for the 2-d hard disk system 
in a 1-d modulating potential showed, in fact, several universal features consistent with the 
predictions of FNR theory. Non universal predictions, namely the phase diagram, on the 
other hand, are difficult to compare because the FNR approach (like KT theory) is expressed 
in terms of the appropriate elastic moduli which are notoriously time-consuming to compute 
near a continuous phase transition. Diverging correlation lengths and times near the phase 
transition point further complicate an accurate evaluation of the non universal predictions of 
the theory. 

In this Letter, we calculate the phase diagram of a 2-d hard disk system with a modulating 
potential (see Fig.^) following a Monte Carlo renormalization approach proposed recently [15] 
for the XY-model. The twin problems of diverging length and time scales are eliminated by 
simulating a constrained system which does not undergo a phase transition! This is achieved by 
rejecting all Monte Carlo moves which tend to distort an unit cell in a way which changes the 
local connectivity [16]. The percentage of moves thus rejected is a measure of the dislocation 
fugacity [16]. This, together with the elastic constants of the dislocation free lattice obtained 
separately, are our inputs (bare values) to the renormalization flow equations [5] to find out the 
melting points and hence the phase diagram. Our results (Fig. ^) clearly show a modulated 
liquid (ML) — > locked floating solid (LFS) — > ML re-entrant transition with increase in the 
amplitude (Vb) of the potential. In general, we find, the predictions of FNR theory to be valid. 
The location of the phase transition as evaluated within this theory with our inputs, show 
excellent agreement with earlier simulations [7] throughout the rj — {3Vq plane (j3 = 1/kgT 
with kg = Boltzmann constant and T — temperature, rj — packing fraction). We discuss our 
calculations in detail below. 
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Fig. 2 - This cartoon shows a typical hard disk system. The dashed lines indicate minima of external 
modulating potential (3V(y) = /3Vo cos(27Tj//do). a x is the lattice parameter and a y indicate the 
average separation between two layers along y-direction perpendicular to a set of close-packed planes. 
For a perfect triangular lattice a y = y3a x /2. The modulating potential is commensurate with the 
lattice such that do = a y . 



The bulk system of hard disks where particles i and j, in 2-d, interact via the potential Vij = 
for \rij\ > d and Vij — oo for |ry-| < d, where d is the hard disk diameter and r»j = rj — 
the relative position vector of the particles, is known to melt [16, 18, 19] from a high density 
triangular lattice to an isotropic liquid with a narrow intervening hexatic phase [16,17,19]. 
Simulations [19], experimental [20] and theoretical [21] studies of hard disks show that for 
rj > .715 the system exists as a triangular lattice which transforms to a liquid below r\ — .706. 
The small intervening region contains a hexatic phase predicted by the KTHNY theory [17] 
of 2-d melting. Apart from being easily accessible to theoretical treatment [22] , experimental 
systems with nearly "hard" interactions viz. stcrically stabilized colloids [20] are available. 
The hard disk free energy is entirely entropic in origin and the only thermodynamically relevant 
variable is the number density p = N/V or the packing fraction rj = (7r/4)pd 2 . 



Theory. - In presence of a periodic external potential, the only other energy scale present 
in the system is the relative potential [23] strength (3V$. A cartoon of the system considered 
for our study is given in Fig. [3] For a solid in presence of a modulating potential f3V(y) (Fig.0) 
displacement mode u y becomes massive, leaving massless u x modes. After integrating out the 
u y modes the free energy of the LFS may be expressed in terms of u x and (3V$ dependent 
elastic moduli [5], namely, the Young's modulus K and shear modulus ji, 



He 
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dy 



(1) 



Similar arguments [5] show that among the three sets of low energy dislocations available 
in the 2-d triangular lattice, only those (type I) with Burger's vector parallel to the line of 
potential minima survive at large (3V$. Dislocations with Burger's vector pointing along the 
other two possible close-packed directions (type II) in the 2-d triangular lattice have larger 
energies because surrounding atoms are forced to ride the crests of the periodic potential [5]. 
Within this set of assumptions, the system therefore shares the same symmetries as the XY 
model. Indeed, a simple rescaling of x — > ^/JIx and y — > V~Ky leads this free energy to 
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the free energy of the XY- model with spin- wave stiffness K xy = y/Kjla f /Air and spin angle 
6 = 2~Ku x /a x . The corresponding theory for phase transitions can then be recast as a KT 
theory [14] and can be described in the framework of a two parameter renormalization flow for 
the spin- wave stiffness K(l) and the fugacity of type I dislocations y'(l), where I is a measure 
of length scale as I = \n(r/a x ), r being the size of the system. The flow equations can be 
expressed in terms of x' — (ixK xy — 2) and y' — At: exp(—j3E c ) where E c is the core energy of 
type 1 dislocations which can be obtained from the dislocation probability [16]. Keeping upto 
next to leading order terms in y' the renormalization group flow equations [15,24] are, 

UA 12 12 I 

■*=-"-« x 

f--*V + f»*. (2) 

Flows in I generated by the above equations starting from initial or "bare" values of x' and 
y' fall in two categories. If, as I — > oo, y' diverges, the thermodynamic phase is disordered 
(i.e. ML), while on the other hand if y 1 vanishes, it is an ordered phase (a LFS) [5]. The two 
kinds of flows are demarcated by the separatrix which marks the phase transition point. For 
the linearized equations the separatrix is simply the straight line y' — x' , whereas for the full 
non-linear equations one needs to calculate this numerically. 

Simulation results and Discussion. We obtain the bare y' and x' from Monte Carlo 
simulations of 43 x 50 = 2150 hard disks and use them as initial values for numerically solving 
the Eqs. J5J). The bare numbers are relatively insensitive to system size since our Monte 
Carlo simulation does not involve a diverging correlation length. This is achieved as follows 
[15, 16]. We monitor individual random moves of a hard disk (after checking for overlaps with 
neighbors) and look for distortions of the neighboring unit cells. If in any of these unit cells 
the length of a next nearest neighbor bond tends to become smaller than a nearest neighbor 
bond, the move is rejected. The probabilities of such bond breaking moves are however 
computed by keeping track of the number of such rejected moves. One has to keep track of 
three possible distortions of the unit rhombus (see Fig. |31 inset) with measured probabilities 
Pmi,i — 1,3. Each of these distortions involves two dislocation- antidislocation pairs which, 
we assume, occur independently. The probabilities for occurrence of the dislocation pairs 
themselves Pdi (Fig. |3J) which are proportional to the square of the fugacities, can then be 
computed easily eg. Pdi = (P m 2 + Pm3 — Pmi)- Finally, we obtain the required dislocation 
fugacity y' by normalizing (Pdi) 1 / 2 to the known [16] dislocation fugacities for (3Vq = at any 
7]. The same restricted Monte Carlo simulation can be used to find out the stress tensor, and 
the elastic moduli from the stress-strain curve, following the method described by Farago et. 
al. [25]. The errors in 'bare' elastic moduli are at worst within a percent. 

In Fig. 0| we have plotted the bare values of x' and y' for r\ — .7029 along with the separatrices 
for the linearized and the non-linear flow equations (Eq. The line of initial conditions is 
seen to cross the non-linear separatrix twice (signifying re-entrant behaviour) while crossing 
the corresponding linearized separatrix only once at high potential strengths. The phase 
diagram (Fig. ^) is obtained by computing the points at which the line of initial conditions 
cut the non-linear separatrix using a simple interpolation scheme. It is interesting to note 
that within a linear theory the KT flow equations fail to predict a RLIF transition. 

Further, comparing with previous computations [7] of the phase diagram for this system 
(also shown in Fig.^) we find that our results agree at all values of r\ and (3Vq. This validates 
both our method and the quantitative predictions of Ref. [5] . The effect of higher order terms 



Debasish Chaudhuri and Surajit Sengupta: A NUMERICAL renormalization GROUP study OF LASER INDUCED FREEZING5 




1000 



PV 

Fig. 3 - Inset shows the unit cell and three possible bond-breaking moves for a triangular solid in 
the orientation shown in Fig. |2] Arrows show the Burger's vectors for the associated dislocations. 
Dislocations with Burger's vector parallel to the horizontal direction are type I dislocations. The 
probabilities for these moves are P m i(top), P m 2 and P m 3- In the main plot the O symbols correspond 
to Pdi, the probability for type I dislocations and the + symbols to (Pd2 + Pd3)f2 the probability 
for type II dislocations obtained from the P m % (see text) for various rj values, arrow denoting the 
direction of increasing r?(= .69, .69395, .696, .7, .7029, .705). These probabilities are plotted against 
the potential strength /3Vo- Note that for /3Vo > 1, the probability for type I dislocations is larger 
than that of type II. 



in determining non-universal quantities has been pointed out earlier [15] for the planar rotor 
model but in the present case their inclusion appears to be crucial. In view of the fact that 
the linear flow equation predicts a solid phase at all but the largest values of (3Vq we require 
at least upto next to leading order corrections in flow to obtain meaningful results. Even 
then, we expect our procedure to break down at high packing fractions in the /3Vb — > limit 
where effects due to the cross-over from a KT to a KTHNY [17] transition at (3Vq = become 
significant. Indeed, as is evident from Fig. [3| for (3Vo < 1 the dislocation probabilities of 
both type I and type II dislocations are similar and our process (which involves only type I 
dislocations) fails to produce melting as (3Vq — > for 77 > .705 — the solid line in Fig. 1 being 
an extrapolation from our results for smaller 77. This fact is also supported by Ref. [7] where 
it was shown that though at /3Vo = 1000 the scaling of susceptibility and order parameter 
cumulants gave good data collapse with values of critical exponents close to FNR predictions, 
at f3Vo — .5, on the other hand, ordinary critical scaling gave better data collapse than the 
KT scaling form, perhaps due to the above mentioned crossover effects. In the asymptotic 
limit of f3Vo — ► 00 the system freezes above 77 ~ .705 which was determined from a separate 
simulation in that limit. This number is very close to the earlier value 77 ~ .71 quoted in 
Ref. [7]. As expected, the freezing density in the (3Vo — > 00 limit is lower than the value 
without the periodic potential. 
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Fig. 4 - The initial values of x' and y' obtained from the elastic moduli and dislocation probability 
at n = .7029 plotted in x' — y' plane. The line connecting the points is a guide to eye. The arrow 
shows the direction of increase in /3Vb(= .01, .04, .1, .4, 1, 4, 10, 40, 100). The dotted line denotes the 
seperatrix (y' — x') incorporating upto the leading order term in KT flow equations. The solid curve 
is the seperatrix when next to leading order terms are included. In I — > oo limit any initial value 
below the seperatrix flows to y' = line whereas that above the seperatrix flows to y — > oo. The 
intersection points of the line of initial values with the seperatrix gives the phase transition points. 
The plot shows a freezing transition at f3Vo = .035 followed by a melting at (3Vo = 38. 

Conclusion. In conclusion, using FNR theory wc calculate the phase diagram for a 
2-d system of hard disks under a commensurate modulating potential to find extremely good 
agreement with earlier simulated results. We show that the re-entrance behavior is built into 
the 'bare' quantities themselves. To obtain the correct phase diagram, however, flow equations 
need to be correct at least upto next to leading order terms in the dislocation fugacity. Our 
results, especially for small potential strengths, is particularly sensitive to these terms. Cross- 
over effects from zero potential KTHNY melting transition are also substantial at small values 
of the potential. 
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